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Numerical  simulations  using  the  lattice  Boltzmann  method  (LBM)  are  developed  to  elucidate  the  dynamic 
behavior  of  condensed  water  and  gas  flow  in  a  polymer  electrolyte  membrane  (PEM)  fuel  cell.  Here,  the 
calculation  process  of  the  LBM  simulation  is  improved  to  extend  the  simulation  to  a  porous  medium  like 
a  gas  diffusion  layer  (GDL),  and  a  stable  and  reliable  simulation  of  two-phase  flow  with  large  density 
differences  in  the  porous  medium  is  established.  It  is  shown  that  dynamic  capillary  fingering  can  be 
simulated  at  low  migration  speeds  of  liquid  water  in  a  modified  GDL,  and  the  LBM  simulation  reported 
here,  which  considers  the  actual  physical  properties  of  the  system,  has  significant  advantages  in  evaluating 
phenomena  affected  by  the  interaction  between  liquid  water  and  air  flows.  Two-phase  flows  with  the 
interaction  of  the  phases  in  the  two-dimensional  simulations  are  demonstrated.  The  simulation  of  water 
behavior  in  a  gas  flow  channel  with  air  flow  and  a  simplified  GDL  shows  that  the  wettability  of  the 
channel  has  a  strong  effect  on  the  two-phase  flow.  The  simulation  of  the  porous  separator  also  indicates 
the  possibility  of  controlling  two-phase  distribution  for  better  oxygen  supply  to  the  catalyst  layer  by 
gradient  wettability  design  of  the  porous  separator. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Water  management  is  essential  to  improve  the  performance 
of  polymer  electrolyte  membrane  (PEM)  fuel  cells.  While  the 
membrane  needs  to  be  fully  hydrated  to  maintain  high  proton  con¬ 
ductivity,  excessive  amounts  of  water  condense  in  the  gas  diffusion 
layers  and  gas  flow  channels  and  prevent  the  supply  of  reactants 
to  the  electrodes  under  high  current  density  conditions.  This  phe¬ 
nomenon  causes  mass  transport  limitations  and  deteriorates  the 
cell  performance.  The  objective  of  this  paper  is  to  evaluate  the 
dynamic  behavior  of  condensed  water  and  gas  flows  in  a  PEM  fuel 
cell  using  numerical  simulations,  and  an  advanced  lattice  Boltz¬ 
mann  method  for  two-phase  flow  with  large  density  differences 
was  developed. 

Understanding  liquid  water  behavior  in  PEM  fuel  cells  is  of 
considerable  practical  significance.  Several  studies  have  been 
conducted  on  two-phase  flows  in  PEM  fuel  cells.  Theoretical  one¬ 
dimensional  models  for  liquid  water  transport  in  gas  diffusion 
layers  (GDLs),  where  liquid  water  is  controlled  by  capillary  forces 
depending  on  the  structure  and  wettability  of  GDL,  have  been 
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reported  [1,2],  and  three-dimensional  simulations  using  two-phase 
models  have  been  developed  [3,4].  For  liquid  water  transport  in  gas 
flow  channels,  experimental  investigations  have  been  conducted 
using  transparent  fuel  cells  and  it  was  demonstrated  that  the  sur¬ 
face  tension  of  water  and  the  wettability  of  the  GDL  and  gas  flow 
channel  play  a  dominant  role  in  the  liquid  water  transport  [5,6]. 
Some  numerical  simulations  of  the  dynamic  behavior  of  condensed 
water,  which  would  be  strongly  affected  by  the  wettability  of  GDL 
and  gas  flow  channel,  have  been  reported.  The  lattice  Boltzmann 
method  (LBM)  was  applied  and  it  was  established  that  LBM  can  be 
a  powerful  tool  to  estimate  two-phase  flow  in  the  gas  flow  channels 
[7,8].  Two-dimensional  simulation  employing  the  volume  of  fluid 
(VOF)  method  were  performed  to  investigate  the  dynamic  behav¬ 
ior  of  a  water  droplet  subjected  to  air  flow  in  the  bulk  of  the  gas 
channel  [9].  Recently,  some  simulations  of  condensed  water  dis¬ 
tributions  in  a  GDL  have  been  developed  using  pore-scale  modes, 
e.g.,  pore-network  modeling,  representing  a  porous  medium  at  the 
microscopic  scale  by  a  lattice  of  wide  pores  connected  by  narrower 
constrictions  termed  throats  [10]  and  a  full  morphology  model  rely¬ 
ing  on  decomposing  digital  images  of  the  GDL  with  pore  radius  as 
the  ordering  parameter  at  a  specified  pressure  during  drainage  [  11  ]. 

This  paper  develops  the  numerical  simulation  using  the  lattice 
Boltzmann  method  (LBM)  to  understand  the  dynamic  behavior  of 
condensed  water  and  gas  flow  in  a  GDL  and  a  gas  flow  channel. 


0378-7753 /$  -  see  front  matter  ©  2009  Elsevier  B.V.  All  rights  reserved. 
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Nomenclature 

c  characteristic  particle  speed  (m  s-1 ) 

C/  restricted  velocities  of  particle  ensembles  (m  s-1 ) 

fi  particle  velocity  distribution  functions  for  the  cal¬ 

culation  of  an  order  parameter 
g  gravitational  acceleration  (ms-2) 

gi  particle  velocity  distribution  functions  for  the  cal¬ 

culation  of  a  predicted  velocity 
H  vertical  length  of  simulation  domain  (m) 

L  characteristic  length  (m) 

p  pressure  (Pa) 

Sh  Strouhal  number 

t  time  (s) 

to  characteristic  time  scale  (s) 

At  time  step  during  which  the  particles  travel  across 

the  lattice  space  (s) 
u  current  velocity  (ms-1) 

u*  predicted  velocity  (m  s-1 ) 

U  characteristic  flow  speed  (m  s-1 ) 

x,  y  position  coordinates  (m) 

Ax  spacing  of  the  cubic  lattice  (m) 

Greek  letters 

Kf  constant  determining  the  width  of  the  interface  of 

two-phases 

Kg  constant  determining  the  strength  of  the  surface 
tension 

fi  viscosity  (Pas) 

§  coordinate  perpendicular  to  the  interface  (m) 

p  density  (kg  m-3) 

Po  reference  density  ( kg  m-3 ) 

o  interface  tension  (N  m-1 ) 

r f,  xg  dimensionless  single  relaxation  time 

0  order  parameter 

0o  reference  order  parameter 

Superscript  and  Subscript 
eq  equilibrium  state 

in  inflow 

G  gas 

L  liquid 

S  solid 

a ,  fi  Cartesian  coordinates 


The  LBM  with  the  simple  algorithm  has  a  number  of  advantages  as 
stated  below,  but  is  attended  with  much  difficulty  in  maintaining 
continuity  at  the  interface  to  simulate  two-phase  flows  with  large 
density  differences  like  condensed  water  and  air  in  a  fuel  cell.  Here, 
the  calculation  process  of  the  simulation  was  improved  to  extend 
the  simulation  in  a  porous  medium  like  a  GDL,  and  stable  and  reli¬ 
able  simulation  of  two-phase  flows  with  large  density  differences 
in  the  porous  medium  was  established.  Using  the  improved  simu¬ 
lation,  the  applicability  of  the  LBM  and  appropriate  conditions  to 
simulate  liquid  water  behavior  in  the  GDL  are  discussed,  and  the 
significant  advantages  of  this  simulation,  which  can  evaluate  actual 
physical  properties,  were  identified.  Further,  two-dimensional  sim¬ 
ulations  demonstrated  examples  of  the  two-phase  flow  affected  by 
the  interaction  between  liquid  water  and  air  flows  in  a  fuel  cell.  The 
effect  of  the  wettability  of  the  gas  flow  channel  on  the  two-phase 
flow  behavior  was  investigated,  and  the  possibility  of  controlling 
the  two-phase  distribution  using  a  porous  separator  without  chan¬ 
nels  is  presented. 


2.  Method  of  simulation 

The  LBM  simulates  mass  and  heat  transport  phenomena  by 
tracking  movements  of  particle  ensembles  where  the  velocities 
are  restricted  by  a  finite  set  of  vectors.  The  particle  population 
is  expressed  by  distribution  functions,  and  the  time  evolution  of 
the  distribution  functions  is  calculated  by  a  simple  law  of  col¬ 
lision  and  transition,  and  it  is  shown  that  macroscopically  the 
LBM  is  equivalent  to  a  continuity  equation  and  the  Navier-Stokes 
equations  for  incompressible  fluids.  Additionally,  introducing  the 
interaction  of  the  particles  in  the  equation  makes  it  possible  to 
simulate  multi-phase  flow.  Because  of  the  simple  algorithm,  the 
LBM  has  a  number  of  advantages:  flexibility  for  complex  boundary 
geometries,  simplicity  for  parallel  computing  and  accurate  mass 
conservation.  In  multi-phase  flow,  tracking  interfaces  is  not  needed 
and  distinct  interfaces  are  maintainable  without  any  artificial 
treatments.  To  simulate  condensed  water  behavior  in  the  three- 
dimensional  gas  flow  channels  of  a  PEM  fuel  cell,  the  advanced 
LBM  proposed  by  Inamuro  et  al.  [12]  was  applied  and  the  authors 
have  confirmed  that  two-phase  flows  with  large  density  differ¬ 
ences,  density  ratios  up  to  1,000,  can  be  calculated  [7,8].  The 
greatest  advantage  of  this  LBM  is  that  it  can  evaluate  the  interac¬ 
tion  between  the  gas  flow  and  condensed  water  in  the  fuel  cell 
where  all  the  properties  and  conditions  such  as  densities,  viscosi¬ 
ties,  an  interface  tension,  wettability  and  flow  velocities  can  be 
simulated. 

In  the  model,  the  non-dimensional  variables  are  defined  by  a 
characteristic  length  L,  a  characteristic  particle  speed  c,  a  charac¬ 
teristic  time  scale  t0  =  L/U,  where  U  is  a  characteristic  flow  speed, 
a  reference  order  parameter  0O,  and  a  reference  density  p0  is  also 
used  [12],  and  “non-dimensional”  is  represented  by  a  circumflex. 
This  paper  uses  a  two-dimensional  9  velocities  model  (2D9V  model) 
and  the  velocities  of  particle  ensembles  are  restricted  to  the  follow¬ 
ing  vectors  c,  (i  =  1,  2, . . .,  9)  in  the  two-dimensional  case  as  shown 
in  Fig.  1  [13]. 

[c1,C2,C3,C4,C5,C6,C7,C8,C9\ 

0  10—1  0  1-1-1  1  1 

“0010-111-1-1  1  J 

Two  particle  velocity  distribution  functions,/;  and  git  are  used. 
Theft  function  is  used  for  the  calculation  of  an  order  parameter  0 
which  distinguishes  two-phases:  0  <  0G  corresponds  to  gas  phase, 
0  >  0£  liquid  phase,  and  0c  <  0  <  0l  the  interface.  The  g,  function 
is  used  for  the  calculation  of  a  predicted  velocity  of  the  two-phase 
fluid  without  a  pressure  gradient.  The  evolution  of  the  particle  dis¬ 
tribution  functions  ft  and  gt  with  velocity  ct  at  point  x  and  time  t 
are  computed  by  the  following  equations. 


Fig.  1.  Lattice  structure  of  two-dimensional  9  velocities  model  (2D9V  model). 
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Here,/jeq  andg^  are  the  equilibrium  distribution  functions,  zy 
and  xg  are  dimensionless  single  relaxation  times,  E;  is  the  associated 
weight  coefficients  presented  below,  Ax  is  the  spacing  of  the  cubic 
lattice,  At  is  the  time  step  during  which  the  particles  travel  across 
the  distance  of  the  lattice  spacing,  p  is  the  density,  jl  is  the  viscosity, 
u  is  the  current  velocity  and  g  is  the  gravitational  acceleration.  The 
subscripts  a  and  (=  x,  y)  represent  Cartesian  coordinates  and 
the  summation  convention  is  used.  The  third  and  last  terms  on  the 
right  hand  side  of  Eq.  (3)  represent  the  effects  of  viscous  stress  and 
gravitation,  respectively. 

The  order  parameter  0  distinguishing  two-phases  and  the  pre¬ 
dicted  velocity  u*  of  the  multi-component  fluids  are  defined  in 
terms  of  the  two  particle  velocity  distribution  functions. 
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i=i 

(4) 

9 

V  ^  ^ 

u  =  l^cigl 

1=1 

(5) 

The  equilibrium  distribution  functions^  andg^  in  Eqs.  (2)  and 
(3)  are  given  by  the  following  equations. 


./P  =H*0  +  Fi 


-  2  d24>  Kf  (  90 
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dka 
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2FKi 

3  'P 
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+  EijGafrp)ciaCW-\FiK^\Vpf 
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where 

e1  =  ±,e2=e3=e4  =  es  =  1,e6=e7=e8=e9  =  ± 


Hi  =  1,H2  =H3  =  ■••  =H9  =  0,  F,  =  -|,Fi  =  3E,(i  =  2,3,  •••,9) 

Kf  is  a  constant  parameter  determining  the  width  of  the  interface 
between  two-phases,  Kg  is  a  constant  parameter  determining  the 
strength  of  the  surface  tension,  and  the  parameters  p0  and  Gap  are 
explained  in  Ref.  [7].  The  interface  tension  d  is  obtained  by  the 
following  equation. 


Here,  f  is  the  coordinate  perpendicular  to  the  interface.  The  first 
and  second  derivatives  are  calculated  using  the  LBM-specific  finite- 
difference  approximations  [12]. 


Because  the  predicted  velocity  it  given  by  Eq.  (5)  is  not  sat¬ 
isfied  by  the  continuity  equation  (V  •  it  =  0),  a  correction  of  it 
is  required.  The  current  velocity  u  which  satisfies  the  continuity 
equation  can  be  obtained  with  the  following  equations. 


Sh 


u-u 

At 


Vp 

P 


(9) 


=  Sh 


V  it 
At 


(10) 


Here,  Sh  =  U/c  is  the  Strouhal  number  and  p  is  the  pressure  of 
the  two-phase  fluid;  note  that  this  definition  gives  the  following 
relationships,  At  =  Sh  Ax,  which  is  represented  by  At  =  A  x/c  with 
dimension  and  means  that  the  particles  travel  across  the  lattice 
space  Ax  during  time  step  At.  This  paper  solved  Eq.  (10)  using  the 
Successive  Over  Relaxation  (SOR)  method.  Details  of  this  model  are 
described  in  a  previous  paper  [7]. 

The  scheme  proposed  by  Seta  and  Takahashi  [14]  was  applied  to 
consider  the  wettability.  In  this  scheme,  the  effect  of  wettability  is 
established  by  the  density  of  the  solid  wall.  Since  the  intermolecular 
force  is  expressed  in  terms  of  density  in  the  LBM,  giving  the  den¬ 
sity  of  a  solid  wall  corresponds  to  giving  the  intermolecular  force 
between  liquid  and  solid  wall.  It  was  confirmed  that  this  scheme 
can  simulate  the  effect  of  wettability  not  only  on  a  flat  surface  but 
also  at  a  corner  inside  a  gas  flow  channel  [8]. 


3.  Developments  to  the  calculation  process 

The  lattice  Boltzmann  method  (LBM)  for  two-phase  flow  with 
large  density  differences  has  been  applied  to  the  simulation  of 
liquid  water  and  air  flow  in  a  PEM  fuel  cell  [7]  and  the  effect  of 
wettability  and  cross-sectional  shape  on  the  liquid  water  behavior 
in  the  gas  flow  channels  was  reported  [8]  by  the  authors.  However, 
there  are  problems  with  the  reliability  of  simulated  results  in  that 
there  is  nonconservation  of  the  mass  of  the  liquid  water.  Further,  dif¬ 
ficulties  to  simulate  the  two-phase  flow  in  a  porous  medium  such 
as  in  a  modified  gas  diffusion  layer  (GDL)  have  been  found.  This 
study  introduced  some  innovations  to  the  calculation  process  and 
their  effectiveness  was  confirmed.  This  paper  presents  an  outline  of 
explanations  of  the  key  improvements  in  the  method  of  solving  the 
Poisson  equation  (10)  and  the  settings  of  the  boundary  conditions 
for  the  pressure  p  on  the  wall  and  the  particle  distribution  function 
fi  at  the  inlet. 

Inamuro  et  al.  solved  the  Poisson  equation  using  an  additional 
LBM  [12],  while  here  the  equation  is  solved  using  the  successive 
over  relaxation  (SOR)  method  as  in  the  case  of  a  conventional 
fluid  simulation  because  of  its  easy  implementation.  The  SOR  is 
a  commonly-used  solution  algorithm  in  conventional  fluid  sim¬ 
ulations  [15].  With  the  SOR  it  is  possible  to  combine  the  LBM 
with  conventional  fluid  simulation,  but  it  is  necessary  to  introduce 
a  suitable  grid  transformation.  In  the  LBM,  co-location  of  a  grid 
where  velocity  and  pressure  are  defined  at  the  same  grid  points, 
as  shown  in  Fig.  2(a)  is  used.  In  the  SOR  method,  the  calculation 
using  the  co-located  grid  induces  spatial  pressure  oscillations.  To 
prevent  the  spatial  oscillation,  the  defined  velocity  grid  points  were 
transformed  as  shown  in  Fig.  2(b),  and  the  SOR  method  used  this 
staggered  grid.  As  the  LBM  treats  the  diagonal  velocities  the  same  as 
the  orthogonal  velocities,  a  formulation  based  on  the  grid  in  Fig.  2(c) 
was  also  conducted. 

The  following  will  present  an  example  of  a  two-dimensional 
simulation  for  the  change  to  the  liquid  water  droplet  in  station¬ 
ary  air  without  gravity,  to  show  the  improvement  that  this  causes. 
The  domain  is  divided  into  60  x  40  square  cells  in  the  x  and  y 
directions  ( Ax  =  2.5  x  10-2),  and  a  liquid  droplet  with  5  Ax  radius 
is  placed  at  the  center.  The  bottom  and  top  of  the  domain  are 
solid  walls,  and  the  left  and  right  sides  are  given  as  free  outflow 
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(a)  Co-locate  grid  (b)  Staggered  grid  1  (c)  Staggered  grid  2 

Fig.  2.  Lattice  configurations  used  in  the  LBM  and  the  SOR  method. 


Fig.  3.  Conventional  and  improved  simulation  results  of  the  mass  change  of  liquid 
water  droplets  in  stationary  air. 


conditions.  The  non-dimensional  parameters  for  LBM  are  iy=  1, 
rg  =  1,  Kf=  0.5(Ax)2,  =  0.092  and  0G  =  0.015,  and  the  proper¬ 

ties  for  the  liquid  water  and  the  air  are  the  same  as  in  the  next 
section.  While  the  liquid  water  droplet  gradually  spread  in  the 
conventional  simulation,  the  droplet  in  the  improved  simulation 
maintained  its  diameter  practically  unchanged.  Fig.  3  shows  a  com¬ 
parison  between  the  conventional  and  the  improved  simulations  of 
mass  change  of  stationary  liquid  water  droplets.  It  is  confirmed  that 
the  mass  conservation  of  the  liquid  water  is  clearly  satisfied  in  the 
improved  simulation. 

In  the  extended  calculation  from  flow  in  a  gas  channel  to  flow 
in  a  modified  GDL,  the  continuity  was  still  broken  and  inappropri¬ 
ate  two-phase  behavior  was  observed  in  some  cases  even  with  the 
improved  method  of  solving  the  Poisson  equation.  This  was  caused 
by  the  boundary  conditions  for  the  pressure  at  the  corner  in  the 
modified  GDL  as  shown  in  Fig.  4(a).  The  boundary  condition  for  the 


►  a< 

1  2 

c 

Solid 

wall 

i 

f  C 

V 

s 

(a)  Pressure  p  at  the  comer 


(b)  Particle  distribution  function  ft 
at  the  inlet 


Fig.  4.  Corrections  when  setting  boundary  conditons. 


pressure  was  set  to  keep  the  pressure  gradient  zero  at  the  wall  in 
solving  the  Poisson  equation,  and  the  pressure  value  on  the  wall 
was  given  by  that  at  the  neighboring  cell  on  the  fluid  side.  Simi¬ 
larly,  the  pressure  at  a  corner  like  point  a  in  Fig.  4(a)  was  given  as 
the  equivalent  of  the  average  pressure  at  the  adjacent  fluid  points  c 
and  d.  This  boundary  condition  sometimes  induced  non-negligible 
errors  in  the  correction  of  the  predicted  velocity  u  using  Eq.  (9). 
This  study  did  not  set  the  pressure  at  the  corner  to  a  certain  value 
but  adjusted  it  in  each  direction  like  ( a-b )  or  (a-c)  to  make  the  pres¬ 
sure  gradients  zero.  This  means  that  point  a  has  different  pressure 
values  to  give  a  zero  pressure  gradient  for  each  direction,  which  is 
a  computational  technique  to  prevent  the  non-negligible  errors.  It 
was  confirmed  that  this  improvement  of  the  boundary  condition 
for  the  pressure  results  in  good  continuity  in  the  modified  GDL. 

Another  problem  of  simulation  in  the  modified  GDL  was  that 
the  liquid  flow  rate  from  inlet  boundary  decreased  when  the  liq¬ 
uid  passed  through  the  GDL.  The  particle  distribution  function  ft 
at  the  first  grid  at  the  inlet  as  shown  in  Fig.  4(b)  was  set  previ¬ 
ous  to  be  the  equilibrium  distribution  function  calculated  by 
Eq.  (6)  using  the  inlet  conditions,  the  order  parameter  0  of  the 
inlet  fluid  and  the  velocity  u.  Flowever,  in  the  case  that  the  pres¬ 
sure  drop  in  the  inlet  fluid  increased,  the  order  parameter  0  at 
the  second  grid  in  Fig.  4(b)  was  increased,  and  the  reverse^  val¬ 
ues  became  larger.  This  induced  the  problem  that  the  mass  flow 
rate  of  the  inlet  fluid  into  the  simulated  domain  decreased.  In  this 
study,  the  inlet  particle  distribution  function  ft  at  the  first  grid  is 
calculated  to  maintain  a  set  flux  0u  between  the  first  and  second 
grids  using  the  values  of  the  second  grid.  This  boundary  condi¬ 
tion  for  the  particle  distribution  function  ft  at  the  inlet  makes  it 
possible  to  control  the  inlet  mass  flow  rate  into  the  modified  GDL 
accurately. 

4.  Results  and  discussion 

4.1.  Liquid  water  behavior  in  a  porous  medium 

Two-phase  relatively  slow  flow  in  a  porous  medium  is  gov¬ 
erned  by  capillary  and  viscous  forces  and  different  flow  regimes 
are  described  by  the  relation  of  these  forces  [16].  The  capillary 
number  Ca  =  upci/a ,  which  represents  the  ratio  of  viscous  to  capil¬ 
lary  forces,  is  an  important  parameter  in  the  behavior  of  condensed 
water  behavior  in  a  GDL:  for  a  typical  fuel  cell  application,  the  capil¬ 
lary  number  Ca  is  of  the  order  of  10-8,  and  two-phase  flow  in  a  GDL 
falls  in  the  regime  of  capillary  fingering  [10].  This  section  discusses 
the  effects  of  capillary  number  Ca  on  the  liquid  water  behavior 
in  a  porous  medium  using  the  LBM  for  the  two-phase  flow  with 
large  density  differences,  and  also  appropriate  conditions  to  sim¬ 
ulate  the  liquid  water  behavior  in  the  porous  medium  of  PEM  fuel 
cells. 
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Fig.  5.  Liquid  water  flow  behavior  in  the  modified  GDL  at  different  inflow  velocities, 
u[n  =  2.0,  0.20,  0.020  ms-1,  Ca  =  2.3  x  10“2,  2.3  x  10“3,  2.3  x  10“4,  H  =  0.1  mm. 


Fig.  5  shows  two-dimensional  simulations  of  liquid  water  flow¬ 
ing  into  a  modified  GDL  with  three  different  water  velocities. 
The  whole  domain  is  divided  into  100  x  50  square  cells  in  the  x 
and  y  directions.  The  vertical  length  of  the  simulated  domain  is 
H  =  0.1  mm.  The  bottom,  y  =  0mm,  corresponds  to  a  micro  porous 
layer  (MPL),  and  there  are  liquid  water  inflows  with  a  velocity  of  u[n 
from  the  20  p,m  wide  pore  at  the  center  of  the  MPL.  The  two  sides, 
x  =  0  and  0.2  mm  are  solid  walls,  and  there  is  a  free  outflow  opening 
at  the  top,  y  =  0.1  mm.  Obstacles  simulating  the  carbon  fibers  in  a 
GDL  are  added,  and  the  porosity  of  the  porous  medium  is  0.8  with  a 
maximum  pore  diameter  of  about  30  p,m.  The  solid  surfaces  are  all 
set  to  hydrophobic  with  the  contact  angle  of  150°,  with  the  order 
parameters  4>s  =  0.025.  The  simulated  GDL  properties,  the  poros¬ 
ity,  the  maximum  pore  diameter  and  the  contact  angle,  are  within 
the  typical  range  of  the  actual  GDLs.  The  GDL  in  this  calculation 
can  be  a  first  stage  model  of  an  example  of  modified  GDLs  in  the 
two-dimensional  simulation.  All  the  properties  of  liquid  water  and 
air  are  the  actual  established  values.  The  density  ratio  of  the  liquid 
water  to  the  air  is  Pl/Pc  =  8 47  (pL  =  997kgm-3,  pG  =  1.18  kg  m-3), 
the  viscosities  of  the  water  and  air  are  iiL  =  8.54  x  10-4Pas  and 
Pc  =  1.86  x  10-5  Pa  s,  and  the  interfacial  tension  between  water  and 
air  is  o  =  7.29  x  10-2  N  m-1 .  The  time  step  At  is  set  to  1.0  x  10-9  s 
and  the  gravitational  acceleration  is  g=0ms-2.  The  other  non- 


dimensional  parameters  for  LBM  are  the  same  as  in  the  previous 
section.  The  finer  time  step  and  grid  size  ensure  better  convergence 
and  accuracy  of  the  simulation  basically,  but  they  induce  longer 
computation  time.  The  appropriate  time  step  and  grid  size  were 
selected  by  analyzing  their  effects  on  the  simulation  of  the  liquid 
water  and  the  gas  flow  in  a  channel  and  a  simplified  GDL.  Fig.  5(a) 
shows  the  simulation  results  with  the  fastest  inlet  velocity  of  liq¬ 
uid  water  u[n  =  2.0  m  s-1 ,  and  the  capillary  number  Ca  =  2.3  x  10-2. 
After  filling  the  first  pore,  the  liquid  water  progresses  through  all 
four  throats  (left  in  Fig.  5(a)).  Then,  it  spreads  out  to  both  sides 
and  grows  into  an  almost  homogeneous  liquid  phase  from  the  bot¬ 
tom  (right  in  Fig.  5(a)).  In  Fig.  5(b),  where  the  inlet  velocity  is 
u[n  =  0.20 ms-1  and  the  capillary  number  Ca  =  2.3  x  10-3,  the  liq¬ 
uid  water  flows  selectively  into  relatively-wide  throats  and  pores 
and  the  direction  of  progress  is  limited  to  two  paths  (Fig.  5(b)). 
With  the  slowest  velocity  u[n  =  0.020  m  s-1  and  the  smallest  capil¬ 
lary  number  Ca  =  2.3  x  10-4  in  Fig.  5(c),  liquid  water  also  progresses 
into  two  throats  in  the  early-stages  (left  in  Fig.  5(c)).  Flowever,  once 
the  largest  path  has  grown  diagonally  upward  left,  the  liquid  water 
in  the  other,  right,  throat  recedes  and  flows  into  the  main  branch: 
finally  only  one  path,  the  one  through  the  largest  throats  and  pores 
is  developed  (right  in  Fig.  5(c)).  This  phenomenon  of  receding  water 
in  a  GDL  is  discussed  in  the  next  section. 

These  simulation  results  show  that  the  selectivity  of  liquid  water 
progression  in  a  GDL  increases  with  decreasing  capillary  number. 
This  is  because  the  effect  of  the  capillary  force  on  water  behavior 
becomes  dominant  over  the  viscous  force.  Although  the  capillary 
number  in  the  simulation  is  much  larger  than  that  in  a  typical  fuel 
cell  application  where  Ca  is  of  order  of  10-8,  the  results  in  Fig.  5(c) 
simulate  the  capillary  fingering  well,  and  is  considered  to  be  close 
to  the  actual  behavior  of  liquid  water  in  a  GDL.  Because  it  requires 
huge  amounts  of  computation  time  to  simulate  the  extremely  slow 
production  of  the  liquid  water  in  the  actual  fuel  cell,  the  following 
simulations  use  0.02  ms-1  as  the  inlet  velocity  of  the  liquid  water 
to  shorten  computation  time. 

Similar  simulations  were  conducted,  with  the  vertical  length 
of  the  simulated  domain  H  changed  from  0.1  mm  to  1.0  mm,  and 
the  time  step  At  1.0  x  10-8  s.  Here,  the  maximum  pore  diameter 
becomes  about  0.3  mm  and  this  simulates  liquid  water  behavior  in 
a  porous  separator  without  channels,  like  that  recently  proposed  as 
an  alternative  to  cells  with  gas  flow  channels  [17,18].  The  results  in 
Figs.  6(a)-(c)  correspond  to  those  at  right  in  Fig.  5.  In  spite  of  the 
10  times  domain  height  H,  the  liquid  water  distributions  are  very 
similar.  These  simulation  results  indicate  that  two-phase  flows  in 
the  porous  medium  with  the  same  capillary  number  are  analogous 
in  the  range  of  relatively-slow  velocities  because  the  length  effects 
on  the  capillary  and  viscous  forces  cancel  out. 

The  LBM  simulation,  considering  capillary,  viscous  and  inertial 
forces,  confirmed  that  it  is  possible  to  simulate  the  similarity  of 
capillary  number  Ca  in  two-phase  flows  in  PEM  fuel  cells.  Further, 
the  LBM  in  this  study,  which  considers  actual  physical  properties, 
offers  significant  advantages  in  evaluating  two-phase  phenomena 
affected  by  the  interaction  between  liquid  water  and  air  flows. 


2.00  x  10'3s  12.0  x  10*3  s  70.0  x  10*3  s 

(a)  ul*  =  2.0  m  s'1  (b)  w/m  =  0.20  m  s'1  (c)  w/"'  =  0.020  m  s'1 

(Ca  =  2.3  x  I  O'2)  (Ca  =  2.3  x  10'3 )  (Ca  =  2.3  x  1 0‘4 ) 


Fig.  6.  Liquid  water  flow  behavior  in  the  porous  separator  at  different  inflow  velocities,  u[n  =  2.0,  0.20,  0.020  m  s_1 ,  Ca  =  2.3  x  10“2,  2.3  x  10“3,  2.3  x  10“4,  H=  1.0  mm. 
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Fig.  7.  Schematic  outline  of  the  computation  domain  for  the  liquid  water  behavior 
in  a  simplified  GDL  with  air  flow  in  the  gas  flow  channel. 


Fig.  10.  Schematic  outline  of  the  computation  domain  for  the  liquid  water  and  air 
flow  behavior  in  a  porous  separator. 


4.2.  Effect  of  air  flow  in  gas  flow  channel  on  the  liquid  water 
behavior 

The  simulation  of  liquid  water  behavior  in  a  GDL  with  air  flow  in 
the  gas  flow  channel  was  conducted  as  an  example  of  a  two-phase 
flow  affected  by  the  interaction  between  liquid  water  and  air  flows 
in  fuel  cells.  Fig.  7  shows  a  schematic  diagram  of  the  computation 
domain:  the  whole  domain  is  divided  into  100  x  35  square  cells  in 
the  x  and  y  directions;  the  vertical  length  of  the  simulation  domain 
is  H=1.0mm;  the  upper  and  the  lower  parts  correspond  to  a  gas 
flow  channel  and  a  simplified  GDL  where  the  channel  height  is  set 
relatively  shallow  to  investigate  effects  of  the  upper  separator  wall 
on  the  liquid  water  because  the  drag  force  of  air  on  liquid  water 
becomes  much  stronger  in  a  two-dimensional  channel  than  in  an 
actual  three-dimensional  channel.  The  uniform  air  flow  and  free 
outflow  condition  are  at  the  left  inlet  side,  x  =  0  mm,  and  at  the  right 
outlet  side,  x  =  2.9  mm,  except  in  the  bottom  parts.  Uniform  flows  of 
liquid  water  are  also  given  at  the  bottom,  y  =  0  mm.  The  velocities  of 
inlet  air  and  liquid  water  are  0.05  m  s-1  and  0.02  m  s-1 .  The  obsta¬ 
cles  simulating  the  simplified  GDL  are  hydrophobic  with  the  contact 
angle  of  150°  (the  order  parameter  0S  =  0.025)  and  the  effect  of 
the  wettability  of  the  separator  is  examined  using  hydrophobic 
(contact  angle  150°,  05  =  0.025)  or  hydrophilic  (contact  angle  40°, 
05  =  0.080)  walls;  At  is  1.43  x  10-7  s,  and  the  properties  of  liquid 
water  and  air,  and  the  other  parameters  for  the  LBM  are  the  same 
as  in  the  calculations  above. 

Fig.  8  shows  the  simulation  results  with  the  hydrophobic  separa¬ 
tor  wall.  The  liquid  water  progresses  selectively  through  the  wider 
pores,  as  in  Fig.  8(a).  When  the  liquid  water  exits  to  the  channel, 
the  liquid  water  body  is  broken  up  and  the  adjacent  water  branch 
in  the  narrow  GDL  pore  recedes,  Fig.  8(a)  and  (b).  This  phenomenon 
has  been  suggested  from  ex  situ  visualization  of  liquid  water  in  a 
GDL  using  fluorescence  microscopy  [19],  and  shows  that  the  LBM 
can  simulate  this  phenomenon.  After  exiting  to  the  channel,  the 


liquid  water  forms  into  a  droplet  and  is  drained  by  the  air  flow  in 
the  channel  (Fig.  8(b)).  Fig.  9  shows  the  results  with  the  hydrophilic 
separator  wall.  The  liquid  water  expelled  from  the  GDL  is  attracted 
by  the  hydrophilic  separator  and  forms  a  water  film  along  the  wall 
(Fig.  9(a)  and  (b)).  This  attracting  of  the  water  produces  a  void  space 
and  air  flow  near  the  interface  between  GDL  and  channel  (Fig.  9(b)). 
This  may  be  expected  to  improve  the  flooding  characteristics. 

The  simulation  results  confirm  that  the  LBM  is  effective  to  sim¬ 
ulate  two-phase  flows  affected  by  the  interaction  between  liquid 
water  and  air  flows,  and  suggests  that  the  wettability  of  the  gas 
flow  channel  is  important  for  the  two-phase  flow  in  a  fuel  cell.  It 
further  suggests  that  control  of  the  wettability  of  the  separator  can 
result  in  a  favorable  distribution  of  liquid  water  and  air  flows. 

4.3.  Liquid  water  and  air  flow  behavior  in  a  porous  separator 
without  channels 

To  ensure  a  favorable  distribution  of  liquid  water  and  reliable  air 
flow,  a  cell  with  a  porous  separator  without  channels  was  simulated. 
The  porous  separator  has  recently  been  proposed  as  an  alternative 
to  cells  with  gas  flow  channels,  and  the  structure  without  chan¬ 
nels  is  expected  to  enable  a  realization  of  uniform  reaction  over 
the  active  area  of  the  membrane  electrode  assembly  [17,18].  Fig.  10 
shows  a  schematic  diagram  of  the  computation  domain:  the  obsta¬ 
cles  simulating  the  porous  separator  and  the  pore  on  the  bottom 
MPL  in  Fig.  6  are  to  the  left,  and  some  further  arrangements  were 
added.  The  whole  domain  is  divided  into  100  x  50  square  cells  in 
the  x  and  y  directions.  The  vertical  length  of  the  simulated  domain 
is  H=  1.0  mm.  The  top  is  a  wall,  and  uniform  air  flow  and  the  free 
outflow  are  provided  at  the  left  inlet  side,  x  =  0  mm,  and  at  the  right 
outlet  side,  x  =  0.2  mm.  A  uniform  flow  of  liquid  water  is  also  pro¬ 
vided  at  the  bottom  pore  of  MPL,  the  velocities  of  inlet  air  and  liquid 
water  are  0.05 ms-1  and  0.02 ms-1.  To  discuss  the  water  and  air 
flow  distributions  in  the  porous  separator,  the  domain  is  divided 
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Fig.  8.  Behavior  of  liquid  water  and  air  flow  in  the  GDL  and  channel  with  the  hydrophobic  separator. 
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Fig.  9.  Behavior  of  liquid  water  and  air  flow  in  the  GDL  and  channel  with  the  hydrophilic  separator. 
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Fig.  11.  Changes  in  water  weight  and  average  air  flow  rates  along  the  separator  in  the  three  regions,  (I)  lower,  (II)  intermediate,  and  (III)  upper  in  the  hydrophobic  porous 
separator. 
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Fig.  12.  Behavior  of  liquid  water  and  air  flow  in  the  hydrophobic  porous  separator,  H=  1.0  mm. 


into  three  regions,  (I)  lower,  (II)  intermediate,  and  (III)  upper,  as 
shown  in  Fig.  10.  Two  cases  with  different  wettabilities  of  the  porous 
separator  were  simulated;  the  solid  surfaces  were  all  hydrophobic 
(contact  angle  150°,  =  0.025)  and  the  upper  shaded  portions 

are  changed  to  hydrophilic  (contact  angle  40°,  4>s  =  0.080);  At  is 


1.0  x  10-7  s,  and  the  properties  of  liquid  water  and  air,  and  other 
parameters  for  the  LBM  are  same  as  the  calculations  above. 

Fig.  11  shows  the  changes  in  water  weight  and  average  air  flow 
rates  along  the  separator  for  each  region  in  the  hydrophobic  porous 
separator.  Before  0.100  s,  the  liquid  water  weight  in  region  I  increase 
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Fig.  13.  Changes  in  water  weight  and  average  air  flow  rates  along  the  separator  in  the  three  regions,  (I)  lower,  (II)  intermediate,  and  (III)  upper  in  the  hydrophobic-hydrophilic 
porous  separator. 


(a)  0. 1 06  s  (b)  0.207  s  (c)  0.240  s 


Fig.  14.  Behavior  of  liquid  water  and  air  flow  in  the  hydrophobic-hydrophilic  porous  separator,  H  =  1.0  mm. 
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first  followed  by  increases  in  region  II,  as  shown  in  Fig.  11(a).  Simul¬ 
taneously,  the  average  air  flow  rates  in  regions  I  and  II  decrease, 
increasing  the  air  flow  rate  in  region  III,  Fig.  11(b),  where  the  ini¬ 
tial  air  flow  rates  at  the  Os  in  each  region  are  very  similar.  Fig.  12 
shows  the  appearances  of  the  areas  of  liquid  water  (black)  and  air 
flow  (uncolored)  in  the  hydrophobic  porous  separator  at  three  time 
points.  Different  from  Fig.  6(c),  the  liquid  water  does  not  move  diag¬ 
onally  upward  left,  but  due  to  the  effect  of  air  flow  it  moves  through 
the  right  throat  (Fig.  12(a)).  Then,  liquid  water  proceeding  upward 
is  broken  up  and  drained  by  the  strong  air  flow  in  region  III,  estab¬ 
lishing  a  dynamic  equilibrium  (Fig.  12(b)  and  (c)).  Thus,  the  increase 
in  the  water  weight  in  region  III  is  suppressed  and  the  oscillations 
in  average  air  flow  rates  in  Fig.  11(b)  after  0.150  s  are  caused  by  the 
breakup  of  the  liquid  water  body. 

Fig.  13  shows  the  changes  in  water  weight  and  average 
air  flow  rates  along  the  separator  for  each  region  in  the 
hydrophobic-hydrophilic  porous  separator.  (The  upper  shaded  por¬ 
tions  in  Fig.  10  are  here  changed  to  hydrophilic.)  Fig.  14  is  the 
appearance  of  the  area  of  liquid  water  (black)  and  air  flow  (uncol¬ 
ored)  at  three  time  points  in  the  hydrophobic-hydrophilic  porous 
separator.  The  liquid  water  weight  grows  in  region  III,  attracted 
to  the  hydrophilic  porous  separator  (Fig.  14(a)).  This  induces  the 
increase  in  water  weight  in  region  III  after  0.106  s  (Fig.  13(a)),  and 
the  average  air  flow  in  region  III  decreases  to  nearly  zero  (Fig.  13(b)), 
forming  a  liquid  water  region  along  the  top  wall  (Fig.  14(b)).  In 
region  II,  the  liquid  water  weight  decreases  and  the  average  air 
flow  rate  increases  drastically  (Fig.  13(a)  and  (b)).  The  liquid  water 
is  broken  up  by  the  strong  air  flow  in  region  II,  and  a  dynamic  equi¬ 
librium  is  maintained  after  0.240  s  (Fig.  14(b)  and  (c)).  It  should  be 
noted  that  the  liquid  phase  along  the  top  wall  is  drained  by  the  air 
flow  in  region  II  continuously  and  maintained  a  certain  amount, 
and  the  dynamic  equilibrium  can  not  be  realized  by  replacing  the 
liquid  water  phase  with  a  solid  wall  in  region  III.  This  air  path  can  be 
anticipated  to  contribute  to  better  cell  performance  because  reac¬ 
tion  gas  is  easily  supplied  to  the  reaction  area  under  the  bottom. 
Overall,  the  simulation  results  show  that  controlling  wettability  of 
the  porous  separator  is  effective  to  realize  an  optimum  two-phase 
distribution,  where  the  air  flow  occurs  between  reaction  area  and 
condensed  water. 

5.  Conclusions 

The  numerical  simulations  using  the  lattice  Boltzmann  method 
(LBM)  for  two-phase  flows  with  large  density  differences  has  been 
developed  to  understand  the  dynamic  behavior  of  condensed  water 
and  gas  flows  in  a  polymer  electrolyte  membrane  (PEM)  fuel  cell. 
The  calculation  process  of  the  LBM  simulation,  e.g.,  the  discretiza¬ 
tion  procedure  in  solving  the  Poisson  equation  and  the  boundary 
conditions  for  pressure  at  corners  and  the  particle  distribution  func¬ 
tion  at  the  inlet,  was  improved  to  extend  the  application  of  the 
simulation  in  a  porous  medium  like  a  gas  diffusion  layer  (GDL). 
With  the  improvements  a  stable  and  reliable  simulation  of  two- 
phase  flows  with  large  density  differences  in  a  porous  medium 


was  possible.  Using  the  improved  simulation,  the  applicability  of 
the  LBM  and  suitably  selected  conditions  to  simulate  liquid  water 
behavior  in  the  GDL  are  discussed,  and  the  significant  advantages 
of  this  simulation,  which  can  consider  actual  physical  properties, 
are  demonstrated.  First,  it  is  shown  that  dynamic  capillary  finger¬ 
ing  can  be  simulated  at  lower  migration  speeds  of  liquid  water  in 
a  modified  GDL,  and  that  the  LBM  simulation,  considering  capil¬ 
lary,  viscous  and  inertial  forces,  simulates  the  similarity  of  capillary 
number,  Ca,  in  two-phase  flow  in  a  PEM  fuel  cell.  Then,  as  examples 
of  the  two-phase  flow  affected  by  the  interaction  between  liquid 
water  and  air  flows,  two-dimensional  simulations  with  air  flow  in  a 
simplified  GDL  with  a  gas  flow  channel  and  in  a  porous  separator  are 
demonstrated.  The  simulation  in  the  simplified  GDL  suggested  that 
the  wettability  of  the  channel  strongly  influences  the  two-phase 
flow  in  a  fuel  cell,  and  the  possibility  that  control  of  the  wettability 
of  the  separator  can  be  used  to  assist  in  an  efficient  distribution  of 
liquid  water  and  air  flows.  The  simulation  in  the  porous  separator 
without  channels  also  showed  that  controlling  the  wettability  of 
the  porous  separator  is  more  effective  to  realize  an  optimum  two- 
phase  distribution  with  the  air  flow  formed  between  the  reaction 
area  and  condensed  water. 
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